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ABSTRACT 

In  this  report  we  show  how  image  formation  or  reconstruction  in  synthetic 
aperture  radar  (SAR)  can  be  viewed  as  the  inversion  of  the  circular  Radon 
transform.  The  advantage  of  viewing  image  formation  in  this  way  is  that 
it  could  be  used  in  situations  where  more  standard  methods  could  fail  such 
as  high  squint  and  ultra-wideband  SAR.  We  examine  previous  work  in  the 
literature  on  circular  Radon  transforms  and  their  inversion.  Next,  we  present 
some  novel  techniques  and  analytic  expressions  for  the  transform  of  some  key 
functions.  We  briefly  consider  motion  compensation.  Finally,  we  propose  a 
number  of  possible  methods  that  could  be  pursued  to  make  new  practical 
image  formation  algorithms. 


20020708  1U 

APPROVED  FOR  PUBLIC  RELEASE 


DEPARTMENT  OF  DEFENCE 


DEFENCE  SCIENCE  I  TECINOLDCY  OICDNISItTION 


fo2- 


DSTO-RR-0211 


Published  by 

DSTO  Electronics  and  Surveillance  Research  Laboratory 
PO  Box  1500 

Edinburgh,  South  Australia,  Australia  5111 

Telephone:  (08)  8259  5555 
Facsimile:  (08)  8259  6567 

@  Commonwealth  of  Australia  2002 
AR  No.  011-855 
August,  2001 


APPROVED  FOR  PUBLIC  RELEASE 


DSTO-RR-0211 


Inverting  the  Circular  Radon  Transform 


EXECUTIVE  SUMMARY 

The  Radon  transform  is  an  integral  along  a  path,  and  in  its  simplest  form  this  path  is 
a  straight  line.  We  will  refer  to  this  form  of  the  transform  as  the  normal  Radon  transform 
(NRT)  because  it  is  usually  written  using  the  normal  equation  for  the  line.  Jakowatz  et  al. 
(based  on  work  by  Munson  et  al.)  showed  the  relationship  between  this  line  integral  and 
spotlight  SAR  image  formation  using  the  polar-to-rectangular  resampling  process.  This 
corresponds  to  the  Fresnel  approximation  to  the  wavefront:  the  radiating  energy  from 
the  radar  is  assumed  to  be  planar.  However,  this  approximation  only  holds  in  restricted 
physical  geometries  and  limits  the  size  and  resolution  of  the  image  that  can  be  formed. 

In  reality,  the  wavefront  radiating  from  the  radar’s  antenna  is  spherical,  with  the  sector 
of  the  sphere  determined  by  the  antenna  beam  pattern  and  geometry.  As  a  consequence, 
the  received  signal  at  a  particular  antenna  position  and  time  since  transmission  is  an 
average  (integral)  of  the  surface  reflectivity  (at  the  radiation’s  frequency)  at  all  points 
where  the  spherical  wavefront  impinged  upon  the  ground  at  half  that  time  (to  account  for 
the  round  trip).  We  will  ignore  the  effect  of  shape  and  duration  of  the  transmitted  pulse 
by  assuming  that  the  received  data  has  been  range  compressed  by  matched  Altering  the 
received  signal  against  the  transmitted  pulse.  The  antenna  position  is  usually  a  function  of 
time  along  a  straight  line,  and  so  this  coordinate  is  referred  to  as  slow  time.  In  contrast,  the 
time  interval  between  transmission  and  reception  of  a  pulse  corresponding  to  a  particular 
range  is  referred  to  as  fast  time. 

If  the  geometry  of  the  imaged  terrain  is  planar,  the  locus  of  its  intersection  with  the 
spherical  wavefront  will  be  the  arc  of  a  circle.  Under  this  assumption,  the  received  radar 
signal  is  an  average  of  the  ground  reflectivity  over  a  circular  arc.  Consequently,  the  radar 
performs  a  circular  Radon  transform  (CRT)  of  the  ground  reflectivity.  Therefore,  forming 
an  image  of  the  ground  reflectivity  in  synthetic  aperture  radar  (SAR)  is  the  process  of 
inverting  the  circular  averages,  or  in  other  words,  inverting  the  circular  Radon  transform  of 
the  received  radar  signal.  This  approach  to  image  formation  provides  an  alternative  view 
to  understanding  the  SAR  image  formation  process  and  can  be  used  to  develop  algorithms 
that  axe  free  of  the  range  curvature  limitations  of  standard  techniques.  The  approach 
used  here  may  also  be  more  intuitive  for  many  readers  unfamiliar  with  radar,  because  the 
concepts  of  doppler  and  phase  are  not  required.  In  addition,  methods  based  on  inverting 
the  CRT  could  be  used  in  situations  where  the  more  standard  techniques  may  fail  because 
the  stationary  phase  approximation  is  inappropriate,  such  as  ultra-wideband  or  foliage 
penetration  SAR,  and  high-squint  SAR  for  tactical  platforms.  Follow-on  work  however 
will  be  needed  to  determine  the  practical  utility  of  the  results  presented  here. 

In  this  report  we  examine  the  existing  body  of  literature  on  inverting  the  CRT  and 
related  methods  in  the  SAR  literature.  In  addition,  we  present  a  number  of  new  methods 
for  inverting  the  CRT  as  the  basis  for  future  research.  We  assume  that  the  radar  is 
monostatic  (co- located  transmitter  and  receiver),  and  that  the  propagation  medium  is 
isotropic,  homogeneous  and  non-dispersive  so  that  the  wave  velocity  is  constant  in  space. 
Then  using  the  Born  approximation,  the  radar  echo  is  assumed  to  be  a  sum  of  single- 
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scattering  objects.  We  also  use  the  “start-stop”  approximation  because  the  sensor  platform 
speed  is  much  less  than  the  wave  velocity.  Furthermore,  we  make  three  assumptions 
that  can  be  removed  at  the  cost  of  complexity.  Firstly,  the  geometry  is  assumed  to 
be  two-dimensional,  i.e.  we  can  ignore  the  height  of  the  radar  above  the  ground  plane. 
Secondly,  the  illumination  of  the  ground  by  the  antenna  in  both  transmission  and  reception 
is  assumed  to  be  uniform.  Finally,  we  assume  an  aspect  independent  ground  reflectivity. 

This  work  was  carried  out  while  the  first  author  was  on  attachment  to  the  Defence 
Evaluation  &  Research  Agency,  Malvern  UK. 
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1  Introduction 

The  Radon  transform  is  an  integral  along  a  path,  and  in  its  simplest  form  this  path  is 
a  straight  line.  We  will  refer  to  this  form  of  the  transform  as  the  normal  Radon  transform 
(NRT)  because  it  is  usually  written  using  the  normal  equation  for  the  line.  Jakowatz  et 
al.  [18]  (based  on  work  by  Munson  et  al.  [23])  showed  the  relationship  between  this  line 
integral  and  spotlight  SAR  image  formation  using  the  polar-to-rectangular  resampling 
process.  This  corresponds  to  the  Fresnel  approximation  to  the  wavefront;  the  radiating 
energy  from  the  radar  is  assumed  to  be  planar.  However,  this  approximation  only  holds 
in  restricted  physical  geometries  and  limits  the  size  and  resolution  of  the  image  that  can 
be  formed  [18]. 

In  reality,  the  wavefront  radiating  from  the  radar’s  antenna  is  spherical,  with  the  sector 
of  the  sphere  determined  by  the  antenna  beam  pattern  and  geometry.  As  a  consequence, 
the  received  signal  at  a  particular  antenna  position  and  time  since  transmission  is  an 
average  (integral)  of  the  surface  reflectivity  (at  the  radiation’s  frequency)  at  ail  points 
where  the  spherical  wavefront  impinged  upon  the  ground  at  half  that  time  (to  account  for 
the  round  trip).  We  will  ignore  the  effect  of  shape  and  duration  of  the  transmitted  pulse 
by  assuming  that  the  received  data  has  been  range  compressed  by  matched  filtering  the 
received  signal  against  the  transmitted  pulse  (see  Section  3.1).  The  antenna  position  is 
usually  a  function  of  time  along  a  straight  line,  and  so  this  coordinate  is  referred  to  as 
slow  time.  In  contrast,  the  time  interval  between  transmission  and  reception  of  a  pulse 
corresponding  to  a  particular  range  is  referred  to  as  fast  time. 

If  the  geometry  of  the  imaged  terrain  is  planar,  the  locus  of  its  intersection  with  the 
spherical  wavefront  will  be  the  arc  of  a  circle.  Under  this  assumption,  the  received  radar 
signal  is  an  average  of  the  ground  reflectivity  over  a  circular  arc.  Consequently,  the  radar 
performs  a  circular  Radon  transform  (CRT)  of  the  ground  reflectivity.  Therefore,  forming 
an  image  of  the  ground  reflectivity  in  synthetic  aperture  radar  (SAR)  is  the  process  of 
inverting  the  circular  averages,  or  in  other  words,  inverting  the  circular  Radon  transform  of 
the  received  radar  signal.  This  approach  to  image  formation  provides  an  alternative  view 
to  understanding  the  SAR  image  formation  process  and  can  be  used  to  develop  algorithms 
that  are  free  of  the  range  curvature  limitations  of  standard  techniques.  The  approach 
used  here  may  also  be  more  intuitive  for  many  readers  unfamiliar  with  radar,  because  the 
concepts  of  doppler  and  phase  are  not  required  [18,  22,  2].  In  addition,  methods  based  on 
inverting  the  CRT  could  be  used  in  situations  where  the  more  standard  techniques  may  fail 
because  the  stationary  phase  approximation  [32]  is  inappropriate,  such  as  ultra-wideband 
or  foliage  penetration  SAR  [16],  and  high-squint  SAR  for  tactical  platforms.  Follow-on 
work  however  will  be  needed  to  determine  the  practical  utility  of  the  results  presented 
here. 

In  this  report  we  examine  the  existing  body  of  literature  on  inverting  the  CRT  and 
related  methods  in  the  SAR  literature.  In  addition,  we  present  a  number  of  new  methods 
for  inverting  the  CRT  as  the  basis  for  future  research.  Like  [36],  we  assume  that  the  radar 
is  monostatic  (co- located  transmitter  and  receiver),  and  that  the  propagation  medium  is 
isotropic,  homogeneous  and  non-dispersive  so  that  the  wave  velocity  is  constant  in  space. 
Then  using  the  Born  approximation,  the  radar  echo  is  assumed  to  be  a  sum  of  single¬ 
scattering  objects.  We  also  use  the  “start-stop”  approximation  because  the  sensor  platform 
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speed  is  much  less  than  the  wave  velocity.  Furthermore,  we  make  three  assumptions 
that  can  be  removed  at  the  cost  of  complexity.  Firstly,  the  geometry  is  assumed  to 
be  two-dimensional,  i.e.  we  can  ignore  the  height  of  the  radar  above  the  ground  plane. 
Secondly,  the  illumination  of  the  ground  by  the  antenna  in  both  transmission  and  reception 
is  assumed  to  be  uniform.  Finally,  we  assume  an  aspect  independent  ground  reflectivity. 

In  the  next  section  we  review  the  basics  of  the  normal  Radon  transform  as  a  starting 
point  for  our  development  of  the  CRT.  Because  the  CRT  is  a  generalisation  of  the  NRT, 
many  of  the  concepts  relevant  to  both  can  be  presented  more  simply  in  the  context  of 
the  NRT  using  well-established  theory.  Having  established  the  basic  Radon  transform 
theory,  we  then  briefly  consider  some  established  SAR  image  reconstruction  methods  in 
Section  3  to  set  the  linkage  with  the  received  radar  signal.  Section  4  states  the  inverse 
CRT  in  concise  terms,  and  Section  5  examines  previous  methods  for  determining  the 
inverse  CRT,  including  the  link  to  the  well-known  range  migration  algorithm  for  SAR 
image  formation.  Section  6  presents  some  new  Fourier  resampling  methods  for  the  ICRT, 
Section  7  derives  properties  of  the  CRT  and  inverse  CRT,  and  Section  8  discusses  various 
options  for  new  practical  methods  for  SAR  image  formation.  Section  9  briefly  discusses 
motion  compensation,  and  conclusions  are  presented  in  Section  10. 


2  Radon  Transform  and  its  Inverse 

2.1  Normal  Radon  Transform 

Typically,  the  Radon  transform  of  the  function  f{x,y),  x,y  €  R.  is  defined  as  a  path 
integral  along  a  straight  line  of  the  function, 

/oo  rOC 

/  f{x,y)6{p-xcos9-ysine)dxdy,  (1) 

-OO  J  — OO 

where  S(-)  is  the  Dirac  delta  function  and  (p,9)  are  the  parameters  of  a  normal  equation 
for  the  line  of  integration.  (See  [28]  for  the  algebra  of  the  Dirac  delta  function.) 

Several  different  approaches  exist  for  inverting  the  Radon  transform.  The  first  one 
we  will  consider  employs  the  Projection  Slice  Theorem  [18],  which  states  that  the  one¬ 
dimensional  Fourier  transform  of  any  projection  pe(p)  =  (7Zf)(p,9)  is  equal  to  the  two- 
dimensional  Fourier  transform  with  respect  to  the  polar  coordinates  of  the  image  f{x,y) 
to  be  reconstructed,  i.e. 

F{ucos9,iysm9)  =  Pg{v).  (2) 

where 

/OO  roc 

/  (3) 

-  OO  J  —OO 

/OO 

-OC 
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Consequently,  by  (2),  we  may  invert  the  Radon  transform  by  use  of  one-dimensional 
and  two-dimensional  Fourier  transforms.  The  standard  Fast  Fourier  Tiansform  (FFT)  on 
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a  uniform  grid  can  be  employed  if  a  polar  to  rectangular  coordinate  transformation  is 
undertaken  [18,  34];  otherwise  non- uniformly  sampled  FFT  techniques  must  be  employed 
[3,  4,  11,  29,  33,  37].  Note  that  the  theorem  can  be  used  in  reverse  order  to  calculate  the 
Radon  transform  of  a  signal  f{x,  y). 

Filtered  backprojection  is  probably  the  most  commonly  used  approach  to  inverting  the 
Radon  transform.  It  is  derived  as  follows.  Starting  with  the  inverse  of  (3), 

/OO  fOO 

/  dk,  dky 

-OO  J  — OO 


and  converting  to  polar  coordinates  we  obtain  [34] 

r*27r  /•OO 


pzir  POO 

/(x,y)=  /  /  uF{ucose,usme)e-^^^''^^^°^^+y^'^^Uude 

Jo  Jo 

noo 

\u\F{i^ cos e,  u  sin Q)e-^^i‘'i^cose+ysin e) 

■OO 

pi:  pOO  /  poo  \ 

Jo  J—oo  \J —OO  / 


(4) 


Typically,  (4)  is  written  as  two  steps,  the  first  of  which  is  a  high-pass  filtering  step  in  the 
frequency  domain. 


peip)  =  I"  1^1 

followed  by  integration  of  the  set  of  projections  y(p,  6)  =  peip), 

pi: 

f{x,y)=  j  g{xcos9  +  ysin9,d)dd 
Jo 

pi:  poo 

—  I  I  g{p,9)5{p  -■  xcos9  —  ysin9)dpd9, 
Jo  J —OO 


(5) 


(6) 


which  is  called  backprojection  and  is  simply  an  integration  along  sinusoidal  paths  of  the 
filtered  projections  g{p^9). 

It  can  be  shown  that  the  backprojection  operation  is  equal  to  half  the  adjoint  Radon 
transform  [9].  To  see  this,  firstly  let  g{p,9)  be  the  normal  Radon  transform  of  some 
function,  then  we  define  the  backprojection  operator  B  by 

pi: 

{Bg){x,y)  =  I  g{xcos9  +  ysin9,9)d9.  (7) 

Jo 

Then  by  the  symmetry  property  g{p,9)  =  g{—p,  -9)  and 

1 

{Bg){x,y)  =  -  J  g{xcos9  +  ysm9,9)d9  (8) 

where  TZ^g  is  the  adjoint  normal  Radon  transform. 
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We  can  see  that  TV  is  the  adjoint  operator  from  the  following  development  [9].  Firstly, 
let  (•,  •)  designate  an  inner  product  in  1R“,  and  in  the  transform  domain  let  (•,  •)p  designate 
an  inner  product  in  R  x  [0, 27r].  Secondly,  let  /(x)  and  /i(x),  x  €  R^,  be  C°°  functions, 
and  g{p,9)  =  72./(x).  Then 


{h,TVg)=  f  /^(x)(7^^6r)(x)o;x 

=  /  h{x)  I  g{x  cos  6  y  sin  6,9)  dO  dx 

Jo 

r  p27T  roo 

=  /  h{x)  /  /  g{p,9)6{p  —  X cos9  —  ysin9)  dpd9 dx 

Jr^  Jo  J —oc 

p2tt  rOO  r 

—  I  /  /  h{x)g{p,9)S{p  —  xcos9  —  ysm9)  dxdpd9 

J 0  J —oo 

p2n  roo 

-  /  (Kh){p,9)g{p,9)dpd9 

Jo  J -oo 


'0  J -oo 
=  {T^h,g)p, 


(9) 


which  confirms  that  7^^  is  the  adjoint  of  TZ  from  the  inner  product  property  of  a  transform 
and  its  adjoint. 

Many  other  schemes  exist  for  inverting  the  normal  Radon  transform.  They  include  the 
“linogram”  and  related  methods  which  employ  a  non-linear  grid  in  the  Radon  domain  to 
reduce  the  two-dimensional  interpolation  of  the  direct  Fourier  slice  theorem  approach  to  a 
one-dimensional  interpolation  in  the  frequency  domain  and  a  one-dimensional  interpolation 
in  the  Radon  domain  [34].  There  are  also  techniques  that  reduce  the  inversion  of  the  Radon 
transform  to  a  problem  in  linear  algebra  through  discretization  [9,  34].  Various  methods 
exist  to  speed  up  the  resulting  algebra  by  using  the  special  structure  of  the  matrices 
involved  [6,  25]. 


2.2  Circular  Radon  Transform 

By  analogy  with  (1),  we  define  the  circular  Radon  transform  of  a  function  to  be  the 
path  integral  of  the  function  along  a  circle  of  radius  t  centred  on  the  point  (u,  0)  on  the 
a:-axis.  This  can  be  written  as 

/OO  roc 

/  /(x,  y)S{t  -  \/{x  -  u)2  -f-  yV  dx  dy.  (10) 

-OO  J  —OO 

The  geometry  of  the  circular  Radon  transform  is  shown  in  Figure  1.  There  are  many  other 
possible  families  of  CRT  depending  on  the  choice  of  centres  for  the  circles:  the  particular 
definition  presented  here  has  been  chosen  because  it  best  corresponds  to  the  standard 
geometry  used  in  SAR. 

From  Papoulis  [28],  the  Dirac  delta  function  S{a{x.y)).  nix^  tj)  —  t  —  \/{x  —  n)^  +  y'^ 
is  a  line  mass  of  density 
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Figure  1:  The  geometry  of  the  circular  Radon  transform. 


Therefore  we  can  rewrite  (10)  as 


g{u,t) 


{Tlcf){u,t)  =  !  j  f{x,y)dxdy. 

J  J  (x  — 


(11) 


The  integral  in  (11)  needs  to  be  carefully  interpreted  as  it  is  an  integral  along  a  line  in 
the  plane  with  respect  to  arc  length.  Let  x  =  u  +  tcosO  and  y  —  tsinO.  Then  from  the 
Jacobian  of  the  transformation  of  variables,  dxdy  =  tdtdO.  Changing  the  intervals  of 
integration  and  eliminating  the  redundant  integral  over  a  single  value  of  t  we  obtain  as 
our  third  equivalent  definition  of  the  CRT, 


r2ir 

g{u,t)  =  {'R.cf){u,t)  =  j  f{u  +  tcos9,tsm9)td9.  (12) 

Jo 

Prom  (10),  the  dual  operator  for  the  CRT  (changing  the  integration  with  respect  to  x 
and  y  to  be  with  respect  to  u  and  t),  called  here  the  circular  backprojection  operator,  is 
defined  to  be 


/OO  pCX) 

/  5(u,  t)6{t  -  v/(a;  -  u)2  +  dt 

-OO  J  —  OO 

/OO 

g{u,yj{x -uf  +  y‘^)du.  (13) 

-OO 

We  will  prove  that  (13)  is  the  adjoint  of  the  CRT  in  Section  2.2.1  below. 

Sometimes  it  is  convenient  to  define  a  semi-circular  Radon  transform  (SCRT)  by 


(7?.c/2/)('U, t)  =  [  f{u  +  tcos9.tsm9)td9. 
Jo 


(14) 
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2.2.1  Adjoint  of  the  Circular  Radon  Transform 

Given  the  CRT  defined  by  (12),  the  dual  operator  defined  by  (13),  and  letting  (•,•) 
and  (-,  •)r  denote  inner  products  in  and  M  x  [0,  oo],  then  we  can  show  that  7Zl  is  the 
adjoint  operator  of  7Zc  as  follows: 

/oo  /-oo  /  r27r  \ 

/  (  /  f  (u  +  t  COS  $,t  sin  0)td0]  g(ii,t)  dt  du, 

■00  Jo  \Jo  ) 

then  setting  x  —  i cos d  and  y  =  tsin 6, 

/OO  poo  poo 

I  /  /(x  -h  u,  y)g{u,  yj y^)  dx  dy  du 

-oo  J  — oo  J  — oo 

/oo  poo  poo 

/  /  /{a;,  y)g{u,  \/{x-  u)2  -f  xf)  dx  dy  du 

/oo  poo  /  poo  _  \ 

/  fix,y)[  giu,  \/ix-  u)2  +  t/)  du  )  dx  dy 

•OO  J -oo  \J -oo  / 

Therefore,  by  the  inner  product  property,  the  dual  operator  TZl  is  the  adjoint  of  TZc- 


3  Reconstruction  Methods  for  Synthetic 

Aperture  Radar 

In  this  section  we  review  simple  methods  for  SAR  image  formation  that  allow  us  to 
establish  the  relationship  between  the  received  radar  signal  and  the  transformed  signal  or 
circular  averages  in  the  CRT  context.  We  follow  Soumekh  [32]  for  the  development  in  the 
first  two  of  the  following  subsections,  with  material  added  to  motivate  the  link  between 
the  radar  signal  and  circular  averages. 


3.1  Time  Domain  Correlation 

Let  p{t)  denote  the  transmitted  radar  pulse.  Then  wdien  the  antenna  moving  along  the 
x-axis  is  located  at  (?i,0),  the  SAR  signature  of  the  ground  point  (a:,,?;,)  is  given  by 

f  2  J {xj  -  u)~  +  tjf 

pit - - - 

I 

where  c  is  the  speed  of  light  in  the  medium.  Without  loss  of  generality,  we  can  assume  that 
fast  time  t  has  been  scaled  so  that  it  represents  the  distance  to  the  point  of  reflectivity, 
effectively  making  ^  =  1.  We  will  use  this  assumption  for  the  remainder  of  this  report. 
Therefore,  the  pulse  becomes 
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In  reconstruction  via  time  domain  correlation,  this  signal  is  correlated  with  the  mea¬ 
sured  SAR  data  in  the  slow-time  fast-time  domain  s{u,t)  to  give  a  measure  of  the  reflec¬ 
tivity  at  the  point  {xi,yi).  Therefore  we  can  write 

=  j  j^s{u,  t)  p*  {xi  -uf  +  yf^  dtdu,  (15) 

where  p*{-)  denotes  the  complex  conjugate  of  p{-)  and  f{xi,yi)  is  the  estimated  ground 
reflectivity.  The  integral  with  respect  to  t  would  normally  be  undertaken  in  the  Fourier 
domain  as  an  integral  with  respect  to  fast  time  frequency  uj  using  Parseval’s  theorem, 

f{xuyi)=  [  [  s{u,  u)  du  du, 

Ju  Joj 

where  P{u})  is  the  Fourier  transform  of  the  transmitted  radar  pulse. 

How  does  s{u,t)  relate  to  the  CRT  of  f{x,y)7  This  becomes  clear  after  the  following 
section  where  we  show  Soumekh’s  backprojection  form  of  SAR  image  reconstruction. 


3.2  Backprojection 

Let  Sm  denote  the  fast-time  matched  filtered  {i.e.  range  compressed)  SAR  signal  given 
by 


Sm{u,t)  =  s{u,t)*p*{-t). 


(16) 


The  underlying  idea  of  fast-time  matched  filtering  is  to  try  and  undo  the  effect  of  convo¬ 
lution  of  the  returns  with  the  radar  pulse  p{t).  Therefore,  in  principle  one  would  do  this 
matched  filtering  with  a  function  q{t)  whose  Fourier  transform  Q{u})  is  a  reasonable  approx¬ 
imation  to  For  the  special  case  where  p{t)  is  an  ideal  chirp  with  infinite  extent,  then 
q(t)  =  p*{-t).  However,  for  other  radar  pulses  p*{—t)  is  not  such  a  good  approximation 
to  the  desired  inverse  filter.  With  the  appropriate  choice  of  q{t),  Sm{u,t)  =  s{u,t)  *q{t) 
would  be  close  to  the  true  unknown  returns  s  that  would  have  been  collected  by  an  ideal 
radar  that  emitted  a  delta  function  pulse. 

Following  (16),  we  can  write 


{xi  -  u)2  -h =  j^s{u,t)p*  -  sJ{xi-uY  +  y‘f^  dt  (17) 


(Xj  -  u)2  -h; 

and  so  we  can  rewrite  the  time  domain  correlation  (15)  as  the  backprojection 


f{xi,yi)  =  J  (u,  {xi  -  uY  +  du. 


(18) 


Consequently,  the  estimated  reflectivity  /  at  the  point  [xi,yi)  is  given  by  the  sum  over 
slow-time  of  all  the  fast-time  points  that  correspond  to  echoes  from  that  point.  This 
corresponds  to  integrating  along  a  hyperbolic  path  in  the  transform  or  slow-time  fast-time 

domain,  because  for  fixed  (xj,  j/,),  the  function  t  =  (x,;  —  u)^  +  yf  describes  a  (conjugate) 

hyperbola  centred  on  (xj,  0)  in  the  (n,  t)  plane  (concave  upwards)  with  asymptotes  of  slope 


±1. 
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3.3  Link  to  the  Circular  Radon  Transform 


The  linkage  between  the  circular  Radon  transform  and  reconstruction  via  backprojec- 
tion  and  hence  time  domain  correlation  of  the  previous  subsections  can  now  be  explained. 
We  can  write  the  received  radar  signal  as  the  integral  over  the  product  of  the  combined 
antenna  pattern,  gain,  range  attenuation,  and  reflectivity  function  with  the  pulse  to  give 
(38], 


s{u,t)  =  j  J  A{u,  X,  y)f{x.  y)p{t-  \/{x-  -0)2  +  dxdy. 


where  A{u,x,y)  is  the  combined  antenna  pattern,  gain  and  range  attenuation  function, 
and  f{x,y)  is  the  reflectivity  function.  Next  we  assume  that  the  reflectivity  function  is 
non-zero  only  in  the  main  antenna  beam  so  that  the  frequency-independent  antenuation 
function  is  unity,  giving 


s{u,t)  =  /  /  f{x,y)p(t-  y {x  -  u)2  -h  yA  dxdy.  (19) 

JyJx  ^  ^ 

Then  substituting  (19)  into  (17)  and  letting  t'  =  \/{x  —  u)"^  +  y^  and  t"  —  yj (x'  -  u)^  -t-  y'^ 
we  obtain 


{uj')=  [  [  [  f{x',y')p{t-t")dx'dy' 

Jt  Uy'  Jx' 

=  [  f  f{x',y')h{u,t',t")dx'dy' 

Jv'  Jx' 


p*  (t  —  t')  dt 


where  h{u,t',t")  is  given  by 


(20) 


h(u,  t")  —  J p[t  —  t")  p*  [t  —  t')  dt 

and  can  be  considered  to  be  the  point  spread  function  of  the  radar.  For  a  chirp  pulse 
h{u,  t\  t")  is  concentrated  at  t' . 

We  can  now  make  the  following  observations.  Firstly,  at  the  radar  sampling  points 
(x,  y),  the  point  spread  function  and  the  delta  function  d{t  -  \/{x-  uY  -I-  y^) 

will  coincide,  so  from  (20)  and  (13)  the  fast-time  matched  filtered  SAR  signal  ,s„,(u,t) 
is  the  function  y(n,  t)  of  circular  averages,  or  circular  Radon  transformed  data,  to  the 
limit  of  the  radar  resolution.  Secondly,  the  backprojection  (18)  is  identical  to  the  circular 
backprojection  (13).  Thirdly,  /(x,  y)  is  not  the  target  or  ground  reflectivity  function  /(.r,  y) 
because  the  adjoint  is  not  equal  to  the  inverse  ,  although  it  can  be  a  reasonable 
approximation  to  it  in  situations  where  the  scatterers  are  isolated  points  such  as  in  sonar. 


4  The  Problem 
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The  SAR  image  reconstruction  problem,  given  the  data  (j{uJ)  consisting  of  integrals 
of  f{x,y)  over  circular  arcs  in  the  plane  centred  on  the  x  axis,  is  to  find  a  closed  form 


DSTO-RR-0211 


solution  for  /.  In  SAR  terms,  g  represents  the  range-compressed  baseband  received  signal 
in  fast-time  t  and  slow-time  u  coordinates  as  discussed  in  the  previous  section.  As  in  (11), 
we  may  write  the  relationship  between  /  and  g  as 


giu,t) 


{ncf){u 


J  J  (x—u 


(x— u)^+J/2zzt2 


f{x,y)dxdy. 


The  geometry  is  shown  in  Figure  1. 

We  shall  interpret  (21)  here  as  the  equivalent  form  shown  in  (12), 


(21) 


f{u  +  tcosd,tsm6)td0. 


(22) 


When  g  is  data  collected  from  a  SAR  system,  then  (22)  should  more  properly  be  [14], 

x27r 

g{u,t)  —  I  h{t)fj,{0)f  {u  +  t  cos  6,  t  sin 6)td0,  (23) 

Jo 

where  h{t)  is  the  range  attenuation  factor  under  the  omnidirectional  isotropic  scattering 
model,  and  /x(0)  represents  the  antenna  diagram  assuming  no  frequency  dependence.  We 
will  assume  that  the  range  attenuation  factor  has  been  compensated  for  in  the  processing 
to  form  g{u,  t)  and  so  can  be  neglected  here.  We  will  also  assume  that  the  target  is  in  the 
main  beam  so  that  the  antenna  diagram  can  also  be  neglected. 

We  may  assume  that  f{x,y)  =  0  for  y  <  0  or  that  f{x,y)  =  f{x,  —y)  as  convenient 
because  either  the  SAR  points  to  one  side  of  the  platform  or  for  low  frequency  SAR,  two 
antennas  are  to  used  to  distinguish  the  left-  and  right-hand  sides  [15].  Note  here  that  the 
CRT  annihilates  functions  fix,y)  that  are  odd  in  y  (i.e.  f{x,y)  =  —f{x,—y)),  and  this 
property  has  important  consequences.  In  particular,  the  inversion  formulae  typically  give 
reconstructions  /  that  are  even,  so  if  we  know  f{x,  y)  =  0  for  y  <  0,  we  will  need  to  multiply 
reconstructions  by  a  factor  of  two  to  give  the  correct  answer  in  the  upper  half  plane. 
Furthermore,  many  of  the  changes  of  variables  to  be  considered  shortly  involve  square 
roots,  and  so  implicitly  restrict  domains  of  definition  or  make  the  function  symmetric. 


4.1  Invertibility  of  the  Circular  Radon  Transform 

The  Radon  transform  integrating  over  all  circles  in  the  plane  is  clearly  overdeter¬ 
mined,  because  the  set  of  all  circles  has  three  dimensions,  which  is  one  greater  than  the 
plane.  The  question  then  is  what  subsets  of  circles  in  the  plane  are  invertible.  Quinto 
[30]  considers  the  case  of  all  translations  of  a  circle  of  fixed  radius  as  well  as  all  circles 
centred  on  a  circle  and  shows  they  are  invertible.  Agranovsky  and  Quinto  [1]  show  that 
the  more  general  case  of  CRTs  along  paths  that  are  not  on  the  zero  sets  of  harmonic 
polynomials  are  invertible.  The  CRT  along  the  lines  considered  in  this  report  is  a  special 
case.  Zalcman  [39]  provides  an  excellent  introduction  to  these  sorts  of  considerations  in 
integral  geometry. 
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5  Previous  Methods  for  Inverting  the  Circular 

Radon  Transform 


5.1  Hankel  Transform  Method 

The  inverse  of  the  CRT  can  be  expressed  analytically  in  terms  of  a  Hankel  transform 
as  follows.  This  was  discovered  independently  by  Fawcett  [12]  and  Andersson  [2],  and 
rediscovered  more  recently  by  Soumekh  [31]  and  Milman  [21,  22]  as  a  statement  of  the  uj-k 
migration  (or  range  migration)  algorithm  [8].  We  will  firstly  derive  the  inverse  in  terms  of 
the  Hankel  transform  as  presented  by  Nilsson  in  [26]. 

Firstly,  define  the  CRT  by 

r2iT 

g{u,t)  =  {TZcf){u,t)  =  /  f{u  +  tcos6,tf,m6)d0.  (24) 

Jo 

Note  that  in  (24),  the  term  t  has  been  neglected  from  the  definition  (12).  It  is  simply 
a  scaling  of  the  radar  signal  g{u,t)  so  it  is  possible  to  neglect  it  here  without  loss  of 
generality.  Next,  take  the  Fourier  transform  with  respect  to  slow  time  u. 


„2iTivu 


g{u^  t)  dll 


00 

00 


=  J  f  {u  +  t  cos  9,  tsmO)  dO^  du 


where  denotes  the  Fourier  transform  with  respect  to  the  first  variable  and  the  second 
variable  is  left  intact.  Changing  the  order  of  integration, 

/OO  r27< 

/  g2;rir(„+i  cos^)g-2,r».f  cosSy  ^  ^  ^ 

-OO  J 0 

i'2-k  roo 

^  ^-2^ivtcos9  e2™'(“+‘^“*^)/(u  +  tcos0,tsin^)dud^ 

^0  J  —  OO 

r’27r 


rZir 

=  /  f  sin  9)  d9. 

Jo 


Now, 


/OC 

-OC 

where  y  =  ts\n9.  Then  changing  the  order  of  integration, 

g^^'^Hv,t)  =  d.9 

/OC  r27T 

f^^'^\lKp)dp  /  e-~^’Hrcose+ps\n9) 

-OC  Jo 


Let 


UCOS0  +  psin(9  =  \/ +  p^  sin(^  —  iji) 
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for  some  ip.  Therefore  using  periodicity  of  the  integrand, 


/•2'k 

Jo 


sm(d—tp) 


dd 


p2Tr 

“/o  ' 


—2nity/v^+p'^  sin  9 


d9 


—  2TrjQ{2nt\/v^  +  p^) 

where  Jo(-)  is  the  Bessel  function  of  the  first  kind  of  order  zero.  Consequently, 

/oo 

(w,  p)  Jo  (27rt\/u2  +  p2)  dp 

■00 


Next,  let 


so  that 


and  so  then 


a  =  y/v^  + 


(25) 


dp  = 


a  da 


—  v'^ 


/OO 

Jo{2Trta) 

-00 


y/  a"^  —  v"^ 


a  da. 


(26) 


We  can  use  the  symmetry  in  the  problem  to  reduce  the  interval  of  integration  in  (26). 
Prom  f{x,-y)  =  /(re,  y),  we  have  f^^'^\v,-p)  =  f^^'^\v,p)  and  from  (25)  we  have 
a  =  v/u2  =  |t;|  when  p  =  Q.  Therefore, 

j|t)|  y/a^  - 

Then  from  the  Hankel  transform  pair  of  order  zero  [5], 

POO 

f^°{q)  =2Tr  f{r)Jo{2nqr)rdr 
Jo 

POO 

f{r)  =  27r  /  f^<^{q)Jo{2iTqr)qdq 
Jo 


we  have  that 


Let 


if  o-  -  lul  >  0 

0  otherwise. 


then 


IpI  =  \/cr2  - 


f^^'^\v,p)  =  Yy(^’^°)(u,  Vp^^T^) 


(27) 


where  is  the  Fourier  transform  of  g  with  respect  to  the  first  variable  u  and  zero-order 

Hankel  transform  with  respect  to  the  second  variable  t  of  g{u,t),  given  by. 


/oo  POO 

-00  J 0 


g{u,  Jo(27rtp)t  dt  du. 


Note  that  f^^'^\v,p)  is  implicitly  even  in  p. 
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5.2  Backprojection  Method 


Andersson  [2]  showed  that  the  backprojection  operator  for  the  CRT  (the  adjoint  ex¬ 
pressed  in  (13))  can  be  used  to  implement  in  (27).  We  give  the  following  develop¬ 

ment  from  Nilsson  [26]. 

Again  we  take  the  backprojection  operator  to  be 


/OO 

g{u,  y/ix-  u)2 

-OO 

then  taking  the  two-dimensional  Fourier  transform  of  TZlg  we  obtain 

/OO  poo  /  poo  _  \ 

/  /  g{u,  ^{x-uf  +  y^)  du  dx  dy 

•OO  j -OO  \j -OO  / 

/OO  /  roo  roc  \ 

^27riuv  /  /  \/(3:  -  u)2  -f 

■OO  V*/  -00  -00  / 


(28) 


du 


/OC  /  TOO  ^00  \ 

^2niuv  /  /  V^x2  -h  y2)g27rwg2.n,x  ^y 

■00  -00  V -00  / 


and  letting  x  =  tcosO^  y  = 


/OO  /  ^OO  ^ZTT  \ 

•OO  vio  io  / 

/OO  ^00  /  r 

/  g(u,t)e2--’  / 

-00  do  Vdo 


,2iTit\/ v^+p'^  s\n{S—il’) 


dO  ]  tdt  du 


/OO  ^OO 

/  g{u,t)jQ(2Ttty/v'^  -h  p 

-OO  d 0 

=  +  P^)- 

Therefore,  from  (27)  and  (29), 

/(^’^)(u,p)  =  Tr\p\iTllg)^^'^Hv,p). 


(29) 


(30) 


5.3  Fast  Backprojection  Method 

Nilsson  [26]  has  developed  a  fast  backprojection  method  for  the  operator  (28).  This 
is  similar  to  the  method  developed  in  [20].  Nilsson’s  method  has  been  applied  to  the 
CARABAS-II  low-frequency  ultra-wideband  VHF  SAR  sensor  [35,  36]. 


5.4  Link  to  the  Range  Migration  Algorithm 

Milman  [22]  has  shown  that  the  range  migration  algorithm  (RMA)  [8,  21,  31]  can  be 
viewed  as  a  close  approximation  to  a  Hankel-Fourier  transform  representation  of  the  inverse 
CRT,  similar  to  the  one  presented  in  Section  5.1.  We  will  now  examine  this  development 
in  detail. 
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Starting  with  the  CRT\ 


g{u,t)  =  jj 


x—u)^+y^=t'^ 


f{x,y)dxdy, 


Milman  noted  that  the  hyperbolic  geometry  in  the  (u,  t)-plane  suggests  the  use  of  hyper¬ 
bolic  functions  to  simplify  the  integrals.  Re-arranging  the  circle  of  integration, 


{u  —  x)^ 

- 9 - 1 - 9  ~  R 

yZ 


we  see  that  it  defines  a  rectangular  conjugate  hyperbola  in  the  (u,  t)-plane  centred  on 
(a;,0),  with  parametric  representation  (ysinh^/;  +  2;,ycosh'^),  -0  e  M.  Therefore,  if  we 
define 


y  sinh  ip  =  x  —  u 

then  using  the  identity  cosh^  ip  -  sinh^  ^  =  1  we  obtain 

x  —  t  tanh  Ip  +  u 
y  =  tsechip. 

The  Jacobian  of  the  transformation  {x,y)  — )■  {t,ip)  is  -tsech^  so  dxdy  =  tsechipdip, 
and  we  can  write 

/OO 

f  .{t  tanh  Ip  +  11,  t  sech  ip)t  sech  ip  dip. 

•OO 


Let 


9o{u,t)  = 


9{u,t) 

t 


Then  taking  the  two-dimensional  Fourier  Transform  of  9o{u,t)  we  have 

f'OO  /•OO  /'OO 


Jf,f) 
9o 


'^\v,k)  =  Iff  f{t  tcinh.  Ip  +  u,t  sech  Ip)  sech  Ip  dip  dudt, 

J  —OO  J  —OO  J  —  OO 


and  letting  t  =  ycosh'^  so  that  dy  —  dt  sech  ip, 


/OO  roo  roo 
-00  J  —OO  J  — 


^27ri{vu-{’ky  cosh  xp) 


f{y  sinh  Ip  -i-  u,  y)  dip  dudy 


00  J  — OO  J  —00 


and  letting  u  —  x-  ysinhV'  so  that  du  —  dx. 


/CO  COO  p 

-00  J  — OO  J  —00 


^27ri{ky  cosh  ip—vy  sinh  ip)  ^27Tivx 


f{x,y)  dip  dxdy. 


(31) 


*  Milman  included  an  term  in  his  definition,  where  A  is  the  wavelength  of  the  radiation.  We  do 

not  need  this  term  because  we  assume  that  the  differing  phase  lags  with  respect  to  frequency  have  already 

been  taken  into  account  by  the  radar  to  measure  p  ^ 
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Next,  we  let  ky  =  (this  is  called  the  Stolt  transformation  [8]).  Also  we  define 


k  V 

cosh(/>  =  — ,  smh(f)——,  (peR, 

Ky  rCy 


and  we  can  see  as  expected  that 


cosh^  4>  -  sinm  (p  = 


-  ^2 


=  1. 


Then, 


cosh(V’ 


(p)  =  cosh  r/'  cosh  (p  -  sinh  tp  sinh  (p 
=  ^  cosh  Ip  -  ^  sinh  ij), 

ky  ky 


and  we  can  rewrite  (31)  as 


{V,ky) 


/oo  roc  roc 

■00  j -00  j -00 


(32) 


The  Hankel  functions  of  the  first  and  second  kind  (also  called  Bessel  functions  of  the 
third  kind)  are  given  by  the  integral  representations  [19] 

1 

J _QQ 
1 

//f)(x)  = - : 

TT*  7-00 


So  we  can  rewrite  (32)  as 

/OO  rOO 

/  Hi^\2nyky)e‘^-‘^'^fix,y)dxdy,  (33) 

-OO  J  —00 

after  ignoring  the  (p  term  because  we  are  integrating  over  an  infinite  domain.  Using  the 
Hankel  functions  as  a  kernel,  the  Hankel  transform  pair  of  order  zero  is  here  defined  to  be 

/OO 

Hq  {27rqr)f{r)r  dr 

■oo 

/oc 

H^^\27rqr)f^^o\q)qdq. 

■oo 


Let  fo{x.,y)  =  f{x-,y)/y.  Then  (33)  can  be  rewritten  as 


9o 


(^Uky)  —  2^/, 


iv,ky), 


(34) 


or  in  other  words,  the  two-dimensional  Fourier  transform  of  go  (after  the  Stolt  transfor¬ 
mation)  is  equal  to  the  Fourier  transform  with  respect  to  x  and  the  Hankel  transform  of 
order  zero  with  respect  to  y  of  fo{x,  y).  Therefore,  we  can  write  that 

y  J -oc  J -oo 
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The  asymptotic  expansion  of  Hq^^z)  for  large  2  is  given  by  [17], 
^0  n!r(-n+i)  ’ 


n=0 


and  taking  only  the  first  term  of  the  expansion  we  obtain 


(36) 


(37) 


(Note  that  this  is  equivalent  to  making  the  stationary  phase  approximation  [22].)  Therefore 
we  can  rewrite  (35)  as 

1  +  i 


1  I  •  /*oo  roo 

f(x,y)  « 

T^\J  2i  J —00  J —00 


(38) 


5.4.1  Shift  to  Scene  Centre 


The  origin  of  the  swath  is  usually  shifted  from  the  antenna  position  to  the  centre  of 
the  swath  to  reduce  the  bandwidth  of  a  signal.  To  achieve  this  let  the  centre  of  the  swath 
be  (a;o,  yo),  and  define  the  new  coordinates  x'  —  x  —  xq  and  y'  =  y  —  yo-  Then  we  rewrite 
(32)  as2 

J  —  00  J  —00  J  —00 

where  f'{x',y')  =  f{x'  +  XQ,y'  +  Xq).  Then  following  through  the  steps  of  the  previous 
subsection,  we  obtain 


1  +  ^ _ y 

7rV2  ^/y'  +  y 


y  J -00  J —00 

(40) 


5.4.2  Statement  of  the  Range  Migration  Algorithm 

We  are  now  in  a  position  to  state  the  RMA  from  (40): 

1.  Compute  the  two-dimensional  Fourier  transform  of  the  scaled  signal  go{u,  t)  to  obtain 
g^f'^\v,k). 

2.  “Match  filter”;  compute  g\^'^\v,k). 

3.  Perform  the  Stolt  transform  ky  —  k"^  —  by  interpolation  and  resampling  of  the 
second  variable  k  to  obtain  y/k^ gQ^'^\v ,  ky) . 

4.  Compute  the  scaled  inverse  two-dimensional  Fourier  transform  (40)  of  the  output  of 
the  previous  step  to  obtain  f'{x',y'). 

Consequently,  we  can  see  that  the  RMA  is  an  approximation  to  a  Hankel-Fourier 
transform  formulation  of  the  inverse  CRT. 

^Note  the  difference  from  (27)  of  [22].  The  presence  of  the  cosh()/»  —  4>)  term  means  that  it  is  incorrect 
to  take  an  term  outside  the  integrals. 
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5.4.3  Validity  of  the  Hankel  Function  Approximation 


How  accurate  is  the  approximation  in  (37)?  This  will  determine  whether  the  novel 
techniques  we  develop  later  in  this  report  have  practical  utility,  because  the  range  migra¬ 
tion  algorithm  is  a  widely  understood  public  domain  algorithm,  which  is  believed  to  cope 
completely  with  range  migration  and  be  applicable  to  high  squint  and  ultra-wideband  SAR 
[8]. 

Prom  Lebedev  [19],  the  residual  of  the  approximation  in  (36)  is  given  by 


\rniz)\  < 


r^(«  +  i)v/2R 

TT  n!(2|z|  sin  J)"’*’! 


[  arg  zj  ^  TT  — 


where  S  is  an  arbitrarily  small  positive  number.  When  n  =  0  (taking  the  first  term  in  the 
expansion  (36)),  r„(2)  simplifies  to 


.(^)l  < 


1 


=  O 


(41) 


8|z|(sin^)  ■ 

Therefore,  the  error  in  the  approximation  (37)  is  of  the  order  of  l/{yky),  but  the  constant 

3 

factor  could  be  extremely  large  due  to  the  presence  of  the  (sin (5)  2  in  (41). 

Figure  2  shows  the  magnitude  of  the  residual 

‘2 


ro{z)  H^^\z) 


(42) 


where  for  our  purposes  ^  =  2nyky  from  (35).  Clearly  the  residual  decays  to  zero  rapidly 
as  expected  from  (41).  We  will  now  examine  the  range  of  values  of  z  =  27ryky  that  are 
likely  to  be  encountered  in  practice. 

The  signal  support  of  the  range  spatial  frequencies  ky  is  determined  from  the  fast-time 
and  slow-time  spatial  frequencies  k  and  u,  respectively,  by  ky  =  k^  —  v'^.  The  fast-time 
spatial  frequency  k  has  units  of  radians/m  and  signal  support  given  by  [8] 


k  e  [Anifc  -  B/2)/c,47r(/e  +  5/2)/c], 


where  fc  is  the  radar  centre  frequency  (Hz)  and  B  is  its  transmitted  bandwidth  (Hz).  For 
the  n-th  target,  its  slow-time  spatial  frequency  v  has  support  [32] 

W  e  fin  =  [2—  COS0n{  —  L),  2—  COS  9n{L)\, 
c  c 

where  0„(u)  is  the  aspect  angle  (measured  off  trajectory)  of  the  7i-th  target  when  the 
radar  is  located  at  (n,0),  L  is  half  the  size  of  the  synthetic  aperture,  and  w  =  27r/  is 
the  transmitted  signal’s  instantaneous  angular  velocity  (radians/s).  For  a  given  fast-time 
frequency  k,  the  slow-time  spatial  frequency  signal  support  of  u  for  the  total  echoed  SAR, 
signal  s{u,t),  denoted  fl^,  is  given  by  the  union  of  that  for  all  targets  in  the  target  anvi. 
For  spotlight  mode  and  a  planar  radar  aperture,  [32] 
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Figure  2:  The  magnitude  of  the  residual:  ro{z). 


where  6c  is  the  average  squint  angle  (measured  off  trajectory),  Xc  is  the  mean  range  of 
the  target  area,  Dx  is  the  length  of  the  antenna,  and  A  =  2nc/LJ  is  the  instantaneous 
wavelength  of  the  radiation. 

Using  the  above  equations  for  the  support  of  fast-time  and  slow-time  spatial  frequen¬ 
cies,  it  is  possible  to  show  when  problems  are  likely  to  occur  with  the  RMA.  In  surveillance 
and  mapping  applications,  the  standoff  to  the  swath  is  usually  measured  in  kilometers.  If 
we  take  the  example  given  in  Table  10.1,  p.  410  of  [8],  the  minimum  slant  range  of  the 
imaged  scene  is  y  ==  712.35  m  for  a  radar  with  a  radar  centre  frequency  of  242.2  MHz. 
The  transmitted  bandwidth  is  B  -  133.5  MHz,  the  half  synthetic  aperture  L  -  380.4 
m,  and  we  assume  that  the  antenna  length  Dj,  =  4  m  for  a  broadside  geometry.  Then 
k  €  [7.3576,12.9496]  radians/m,  and  u  €  0*  «  [—5.94043,5.94043]  at  the  minimum  instan¬ 
taneous  frequency.  Therefore  z  =  2'Ky\/k^  —  —  19430.3  in  this  example.  However,  by 

shortening  the  antenna  length  Dx  it  is  possible  to  make  ky  arbitrarily  close  to  zero  which 
puts  one  in  the  regime  where  the  approximation  residual  is  large.  Further  investigation  is 
required  to  understand  when  this  is  likely  to  occur  in  practical  systems,  but  this  simple 
example  shows  that  it  is  possible  in  practice. 


5.5  Deconvolution  Method 

Norton  [27]  in  the  context  of  acoustic  imaging  has  developed  a  method  to  reconstruct 
a  reflectivity  field  from  line  integrals  over  circular  paths.  Following  [27],  we  may  rewrite 
(21)  in  terms  of  Dirac  delta  functions  so  that 

/oo  poo 

/  f{x,y)S 

-OO  J  — oo 


(t  —  dxdy. 


(43) 
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Norton  then  uses  the  property  of  the  Dirac  delta  function  that 

S(z  -  a)  =  2z5{z^  -  a^) 


(44) 


for  3;  >  0  and  a  >  0  to  write 

/CX)  poo 

/  f{x,y)5{f  -  {x -uf -y^)  dxdy.  (45) 

-cc  J  —  00 

This  manipulation  is  not  correct,  however,  because  the  substitution  of  (44)  into  (43)  has 
used  z  =  t  and  a  =  ,  which  is  only  appropriate  if  the  integral  (43)  is  with 

respect  to  t,  but  it  isn’t.  Consequently,  the  equations  Norton  derives  to  invert  (45)  are 
not  the  inverse  of  (43). 

6  Resampled  Fourier  Methods  for  Inverting  the 
Circular  Radon  Transform 

We  now  present  a  novel  method  for  inverting  the  circular  Radon  transform  via  Fourier 
transforms  and  appropriate  resamplings.  We  derive  the  result  by  two  different  methods 
and  show  that  they  are  equivalent;  firstly  directly  via  change  of  variables,  and  secondly 
via  the  projection  slice  theorem. 


6.1  Square  Fourier  Transform  Method 

Firstly,  let  us  define  the  square  Fourier  transform  (FT)  of  the  received  signal  g  by 


G{u,u) 


-i: 


t)  dt, 


(46) 


then  from  (22) 


/oo  r2iT 

/  +  tcos6,tsm9)tdtd9. 

-00  J  0 


Making  the  change  of  coordinates  from  polar  {t,  9)  to  Cartesian  (.?:,  y)  (dx  dy  —  tdf  d9  from 
the  Jacobian)  and  then  translating  in  x  so  that  {u.  0)  is  now  the  origin  gives 


TOO  roo  , 

G{u,u)^2  /  ^f(x,y)dxdy 

J -OO  J -00 

aoc 

-OO 


,2iTiuJu^ 


(47) 


Next,  make  the  change  of  variable  r  =  x^  +  y^  so  that 


X  —  X 

y  =  vr  —  .T^ . 


(48) 
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The  Jacobian  of  this  transformation  is 


d{x,y)  _ 

dx 

dx 

dx 

dr 

i 

1 

0 

1 

1 

d{x,  r) 

dx 

dr 

v/j — x'^ 

2\/r—x‘^ 

to 

1 

to 

so  then 

dx  dy  =  —  ^  dx  dr. 
2vr  —  x^ 

We  define  the  function  k{x,  r)  by 


k{x,r)  = 


f{x,Vr—x‘^) 

Vr—x'^ 

0 


0  <  \x'^  <  r 

otherwise, 


(49) 


where  we  have  assumed  that  f{x,  y)  =  f{x,  -y)  so  we  need  only  consider  the  positive  root 
of  vr  — 


Then  we  can  rewrite  G{u,u)  in  (47)  as 

roo  roo 


G{u,ui)  =  2e 


2iTiu)u'^ 


where 


=  e 


2niu}u^ 


/oo  po 

•00  J -oo 

/oo  roo 

■00  J  —I 


f{x,\/r  -  x^) 
2\/r  —  x^ 


dxdr 


J2wi{—2ijju)x  2-Kiujr  i 


k{x,r)dxdr 


=  K{-2oju,u) 


/OO 

/ 

-oo  J  —  c 


e^^i<^^e^^^^^k{x,r)dxdr 


^  O' 


(50) 


(51) 


Calculating  the  Inverse 


Consequently,  given  g{u,  t),  the  inversion  of  the  CRT  can  be  carried  out  by  the  following 
series  of  resamplings  and  Fourier  transforms.  Firstly,  let  h{u,  r)  be  a  scaled  resampling  of 
g{u,  t)  defined  by 


h{u,T)  = 


9{u,  y/r) 
2^ 


(52) 


This  resampling  converts  the  square  FT  of  g  in  (46)  into  the  standard  FT  of  h,  so  that 


G{u,uj)  =  j 

Next,  we  resample  G  according  to 


„2wiujT 


h{u.,  t)  dr. 


(53) 


(54) 
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and  then  compute  k{x,  r)  by  the  two-dimensional  FT  of  K  by 

/OO  rOC 

/  (55) 

-OO  J  — OO 

which  we  resample  to  give  f{x,y)  by 

f{x,y)=^yk{x,x‘^+i/).  (56) 

Note  that  we  have  used  this  method  to  compute  some  analytic  inverse  CRTs  in  Section  7.3. 


6.2  Reduction  to  a  Normal  Radon  Transform 


Starting  with  (22)  we  make  the  change  of  variable  x  =  u  +  t  cos  B  so  that 


dx  =  ~  t  sin  BdB 


and  using  f{x,y)  =  f{x,  -y),  we  can  rewrite  the  integral  in  x  as 


g{u,t)  = -2  f  f{x,  -{x-  u)^)-^=|:.  .dx. 

Ju+t  -  [x  —  uY 

Now  y  =  sjY  —  {x  —  uY,  and  making  the  change  of  variables  (x,  y)  i — >  (x,  r)  as  in  (48), 
i.e.  r  =  x^  +  we  have  that  \/r  -  x^  =  y/Y  -  (x  -  uY  and  r  =  Y  -u^  +  2xu.  Therefore, 
from  the  definition  of  k(x,r)  in  (49),  we  obtain  that 


2xu)t  dx. 


(57) 


We  recognise  this  as  the  integral  along  the  straight  line  xcos0  -I-  rsin^  =  p  in  the  (x,r) 
plane  where 


cos  9  = 


2u 


Vi-  +  4u^ 


sm 


B  = 


Vi-  +  4^2  ’ 


Y  — 

\/l  -f  Av? 


(58) 


However,  the  integral  is  not  with  respect  to  arc  length  along  the  line:  to  achieve  this  we 
note  that 


dl  =  \/l  -t-  4rY  dx. 


(59) 


which  takes  into  account  the  length  of  the  line,  so  that  the  normal  Radon  transform 
(77.fc)(p,  0)  of  k{x.r)  is  given  by 


{'R.k){p,9)  =  /  k{x,r)d.l 

J (j-.r):  X  oos  S-f-r  sin  S=p 


r  oc 

==  \/l  +  4u“  /  k{x.  /“  —  ?r  -f  2xa)  dx. 
J  -OO 


(60) 
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We  can  extend  the  limits  of  integration  of  (57)  to  (-00, 00)  by  noting  that  k{x,r)  =  0  for 
I  a;  I  >  r,  r  <  0.  Therefore,  we  can  rewrite  (60)  as 

(61) 

Nilsson  [26]  has  developed  a  very  similar  expression  to  this  to  relate  the  CRT  to  backpro- 
jection  along  straight  lines. 


By  appealing  to  the  projection  slice  theorem,  we  can  write 


/oo 

e^'^‘^^P{'Rk){p,0)dp 

•OO 


-L 


Prom  (58)  we  can  rewrite  this  as 


/oo 

•00 


°°  2^0.  y/l  +  4^2  ^  ^ 

e  v'i+4n2 - - - g{u,t) 


2t 


2t 


\/l  +  A.v? 


dt 


=  e  \/l+4u2  /  g  \/l+4u2  g[u,t)dt 

J  —  OO 


=  e  4wsine  C(u,  wsin^).  (62) 

By  letting  a  =  uj  cos  0  and  ^  =  u  sin  6  we  obtain 

K{a,0)^e-^^WG{-~,p)  (63) 

provided  that  w  ^  0.  Note  that  in  (63)  we  have  obtained  the  same  expression  as  in  (51), 
confirming  the  earlier  result. 


Calculating  the  Inverse  by  the  Normal  Radon  Transform 


Examining  the  derivation  of  the  CRT  and  its  relationship  with  the  normal  Radon 
transform  above  leads  us  to  the  following  series  of  resamplings  and  Fourier  transforms  as 
an  alternative  to  the  ones  in  Section  6.1.  Firstly,  scale  g{u,t)  to  determine  the  normal 
Radon  transform  of  k{p,9)  using  (58)  by  defining  the  resampled  function  h{p,9), 


hip,  9)^ 


-1 


2sin0 


\/  sinS 


+ 


1 

4  tan2  6 


+ 


2  tan  ^  v  sin  0  4  tan^  9 


so  then 


(64) 


ink){p,9)=^h{p,9) 

by  (61).  Next,  we  invert  the  normal  Radon  transform  {7lk){p,9)  to  obtain  k{x,r)  by  one 
of  the  methods  outlined  in  Section  2.1,  i.e.  , 

kix,r)  —  {K~^{Tlk))ix,r). 

Then  we  can  resample  k{x,r)  to  obtain  fix,y)  by 

fix,y)  =  ykix,x‘^  +  y^). 
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7  Properties  of  the  Circular  Radon  Transform 

In  this  section  we  examine  a  number  of  properties  of  the  CRT.  These  include  the  effect 
of  translation  of  the  target  function  in  the  transform  domain,  and  the  CRT  and  inverse 
CRT  (ICRT)  of  some  simple  functions.  These  properties  are  important  because  we  will 
use  them  later  to  suggest  further  new  methods  for  SAR  image  formation. 


7.1  Translations 

Firstly,  note  that  the  translation  in  x  corresponds  to  a  translation  in  u  in  the  transform 
domain  from  the  definition  (22)  so  that 


(7^c/(a:  -  X,:,  y)){u,  t)  =  g{u  -  x^t).  (65) 

where  {TZcf{x,y)){u,t)  =  g{u,t).  Note  that  this  applies  for  all  f{x,y). 

Next,  let  us  consider  the  effect  of  translation  perpendicular  to  the  imaging  trajectory. 
Let  us  consider  a  disc  as  our  target  function  /(x,y),  but  this  time  displaced  from  the 
origin  along  the  y  axis  (Figure  3), 


1  if  +  (y  -  y;)2  ^ 
0  otherwise. 


This  target  function  is  radially  symmetric,  which  simplifies  considerations  for  the  moment. 

Figure  3  also  shows  the  circular  integration  path  of  radius  t  centred  at  (u,  0)  intersecting 
the  disc  centred  at  (0,  y,).  It  is  clear  from  this  diagram  that  the  resulting  path  integral 
is  equivalent  to  that  obtained  from  a  disc  centred  at  the  origin  integrated  along  a  path 
centred  on  (^u^  -|.  y?,0).  Consequently,  for  radially  symmetric  target  functions,  we  may 
write 

{'^cfix.y  -  yi)){uA)  =g{^u'^  +  yf,i)-  (66) 


For  target  functions  that  are  not  radially  symmetric,  clearly  the  target  function  imaged 
at  (0,  yj)  has  undergone  a  rotation  with  respect  to  that  imaged  at  the  origin.  The  rotation 
is  through  an  angle  of  if)  =  —  tan""^  ^  which  changes  with  slow-time  u.  Then  we  can  write 

(7^c/(x,y  -  iyi))i^ij,)  =  y_tan->  +  (67) 

where  g^,{u,t)  is  a  function  of  iJj  as  well  as  u,t,  and  is  the  transform  of  the  rotated  target 
function  given  by 

gji,{u,t)  =  {TZef  {x  cos -ij)  +  y  sin  iJk  -x  sin  iji  +  y  cos  ij))){u,t).  (68) 

Consequently,  to  determine  the  CRT  of  a  target  function  translated  along  the  y-axis,  w(' 
need  to  know  its  CRT  at  all  rotations  V'  when  it  lies  on  the  x-axis. 
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Figure  3:  The  geometry  of  a  disc  displaced  from  the  origin  along  the  y-axis,  perpendicular 
to  the  imaging  trajectory. 


7.2  Analytic  Forward  Transforms 


In  this  section  we  explore  a  number  of  different  analytic  expressions  for  the  forward 
circular  Radon  transform. 


7.2.1  Delta  Function 

Let  our  target  function  be  the  delta  function  located  in  the  {x,y)  plane  at  the  point 
(0,yi).  Then  we  may  write 

f{x,y)  =  S{x)S{y  -  yi), 

where  5(-)  is  the  Dirac  delta  function.  The  CRT  of  the  delta  function  is  then  given  by 


/»z7r 

g{u,t)  =  I  f{u  +  tcos9,tsm9)td6 

Jo 

=  /  5{u  +  tcos9)6{tsm9 -yi)td9, 

Jo 
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u 


Figure  J^:  The  trace  of  the  circular  Radon  transform  of  a  Dirac  delta  function  at  (0,2). 
The  weighting  function  on  the  implicit  curve  in  (70)  is  not  shown. 


then  from  Papoulis  [28],  S{a{t))  =  where  U  are  the  zeros  of  a{t), 


6{e  -  cos-^  =f)  6{e  -  sin-^  f] 
0  t  sin  cos“i  t  cos  sin“^  ^ 

p27r 


t  dO 


f  5 {9  -  cos  ^ —  sin  ^  ] 
Jo  i  t  ‘ 


de 


then  solving  cos 


-1  —u 


sin  ^  for  t  and  taking  the  positive  root  we  obtain 


g{u,t) 


0  otherwise. 


ift  =  ^u2  +  y? 


(69) 


(70) 


Figure  4  shows  a  plot  of  the  trace  swept  out  by  the  CRT  of  the  delta  function.  Note  the 
hyperbolic  shape,  which  is  characteristic  of  the  CRT  of  finite  or  concentrated  objects  in 
the  (x,  y)  plane. 


7.2.2  Disc 

Let  us  consider  a  disc  centred  on  the  origin,  with  its  circumference  given  by 

+  y2  = 
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We  wish  to  determine  the  CRT  of  this  disc  for  irradiating  circles  of  radius  t  centred  on 
the  x-axis  at  the  point  (u,0).  To  do  this,  we  need  to  determine  the  length  of  the  arc  of 
the  circle  that  intersects  the  disc. 

Let  g{u,  t)  denote  the  resulting  CRT  of  the  disc.  In  the  case  where  |u|  >  |r+t|,  the  circle 
will  not  intersect  the  disc,  and  g{u,  t)  =  0  in  such  circumstances  (Figure  5(a)).  When  |ul  = 
\r  + 1|,  there  is  a  non-enveloping  osculation  between  the  disc  and  the  circle  (Figure  5(b)), 
but  because  the  intersection  between  the  two  occurs  at  an  infinitely  small  point,  g{u,  t)  =  0. 
When  |r  -  t|  <  |u|  <  |r  -F  <|,  intersection  occurs  between  the  circle  and  disc  (Figure  5(c)), 
and  the  length  of  the  intersecting  arc  is  given  by  g{u,t)  =  2tcos“^ 
enveloping  osculation  occurs  when  |u|  =  |r  —  t|  and  t  ^  r  (Figure  5(d)),  and  the  entire 
circle  lies  within  the  disc,  i.e.  ,  g{u,t)  =  2iTt.  If  |u|  =  |r  —  t|  and  t  =  r  then  the  circle 
lies  on  the  circumference  of  the  disc  and  g{u,t)  =  2nt.  If  |u|  <  [r  -  tj  and  t  ^  r  then 
the  circle  is  completely  enveloped  by  the  disc  and  g{u,t)  =  2-Kt  (Figure  5(e)).  However,  if 
|u|  ^  1^  -  and  t  >  r  then  the  circle  envelopes  the  disc  and  g(u,t)  —  0,  intersecting  at 
most  at  a  single  point  when  |u|  =  |r  —  t|. 

In  summary,  the  CRT  of  the  disc  centred  at  the  origin  is  given  by 


9{u,t)  =  { 


2(cos-‘  (s^^) 

27rt 


0 


if  |r  -  t|  <  |u|  <  |r  -j- 1\, 
if  |ix|  ^  |r  —  t|  and  t  ^  r, 
otherwise. 


(71) 


Figure  6  shows  a  density  plot  of  the  CRT  of  a  disc  of  radius  one.  Figure  7  shows  the  effect 
of  translation  along  the  y-axis  on  the  transform  of  a  disc  of  radius  one. 


7.2.3  Gaussian  Radial  Basis  Function 

Let  the  Gaussian  radial  basis  function  (RBF)  be  defined  by 

/(x,y)  =  e-"("^+^').  (72) 

Then  the  CRT  of  f(x,y)  from  (22)  is  given  by 

r2TX 

g{u,t)  =  {'R.cf){u,t)  =  j  f{u  +  tcos6,tsmd)td9  (73) 

Jo 

=  27rte-’^(^'+“'^/o(27rtu),  (74) 

where  /o(-)  is  the  modified  Bessel  function  of  the  first  kind  of  order  0.  Figure  8  shows  a 
density  plot  of  the  CRT  of  the  Gaussian  radial  basis  function. 

We  can  check  the  effect  of  translation  of  the  Gaussian  RBF  along  the  y-axis  to  con¬ 
firm  the  relationship  in  (66)  for  radially  symmetric  target  functions.  Figure  3  shows  the 
geometry  of  the  problem.  If  we  denote  the  -  tan“^  ^  rotated  coordinate  system  (x',y') 
centred  on  (0,  y,;)  in  the  (x,  y)  plane  with  x'  axis  passing  through  the  x-axis  at  x  =  u,  then 
we  can  rewrite  the  Gaussian  RBF  function  as  /(x',  y')  =  K  The  centre  of  our 

integrating  circle  is  located  at  {Ju^  +  yf,0)  in  the  rotated  coordinate  system,  so  then  the 
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CRT  of  f{x,  y)  from  (22)  is  given  by 

g{u,t)  =  J  +  yf  +  t  cos  6,  t  sin 6)tdd, 

which  integrated  gives  that 

g{u,t)  =  {TZcfix^y-  yi)){t,u)  =  27rte-"(''+“'+^?)/o(27rty^u2  +  yf). 

Figure  9  shows  the  effect  of  translation  along  the  y-axis  on  the  transform  of  the  Gaus¬ 
sian  RBF. 


7.2.4  Vertical  Line 

Let  our  target  function  be  the  vertical  line  at  j:  =  0.  Then  we  may  write 

f{x,y)  = 

The  CRT  of  the  vertical  line  is  then  given  by 

/OO  rOO 

/  f{x,y)S{t  -  \/{x  -  u)2  -I-  y^)dxdy 

-OO  J  —  OO 

/OO  roc 

I  6{x)6{t  —  \/{x  -  -I-  y2)  dx  dy. 

-OO  J -OO 

Following  Papoulis  [28],  the  angle  of  the  intersection  of  the  two  line  masses  S{x)  and 
6{t  —  \/{x  —  +  y2)  is  given  by 


sin0  =  \  1  — 


and  so  the  density  at  each  intersection  is  1/ \Jl  —  Therefore, 


if  |u|  <  t 


g{u,t)  =  I  \A^ 

I  0  otherwise. 


7.2.5  Vertical  Line  Segment 


Let  our  target  function  be  the  vertical  line  segment  at  x  =  0  in  the  interval  y  €  [yo,  yi] 
where  0  ^  yo  <  yi-  Then  we  may  write 


/(a.%y) 


S{x)  ifyo<y<yi 
0  otherwise. 


The  CRT  of  the  vertical  line  segment  is  then  given  by  the  expression 


giuj.) 


-7-- if  \/u2  +  yo  ^  t  ^  0/2  -f-  yf 

Vi-75- 

0  otherwise. 


using  the  density  at  the  intersection  calculated  in  Section  7.2.4.  Figure  10  shows  a  density 
plot  of  the  CRT  of  a  vertical  line  segment. 
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(a)  Non-enveloping  non-intersecting:  In]  >  [r-f-  (b)  Non-enveloping  osculation:  lii|  — |r-|-t| 

*1 


(c)  Intersecting:  [r  -  t|  <  |u|  <  |r-i-t|  (d)  Enveloping  osculation:  |u|  =  k  -  t|  and 

t  <  r 


(e)  Enveloping:  |u|  <  |r  - 1|  and  t 


Figure  5:  The  various  cases  under  which  a  circle  can  intersect  a  disc. 


27 


DSTO-RR-0211 


Figure  6:  The  circular  Radon  transform  of  a  disc  of  radius  one  centred  at  the  origin. 


Figure  1:  The  CRT  of  a  disc  of  radius  one  displaced  from  the  origin  along  the  y  axis  to 
be  centred  on  y  =  \/5. 
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7.2.6  Horizontal  Line 

Let  our  target  function  be  the  horizontal  line  at  y  =  yi.  Then  we  may  write 

fix,y)  =  6{y-yi). 


The  CRT  of  the  horizontal  line  is  then  given  by 

/oo  poo 

/  /(a;,y)<5(t  -  y/(x  -  up  -hy^jdxdy 

-OO  j  —  OO 

/oo  poo 

/  Hy  -  yi)Hi  -  +y‘^)dxdy. 

■OO  j  —  00 

Following  Papoulis  [28],  the  angle  of  the  intersection  of  the  two  line  masses  5{y  —  yi)  and 
6{t  ~  \/{x  —  uY  +  is  given  by 

sin(^  -  9)  =  cos6  —  — , 

2  i 

and  so  the  density  at  each  intersection  is  A.  Therefore, 


/  /  2|^  if  t  >  yi 

0^  otherwise. 


7.2.7  Horizontal  Line  Segment 

Let  our  target  function  be  the  horizontal  line  segment  at  y  =  yi  in  the  interval  x  € 
[a;o,a;i]  where  xq  <  xi.  Then  we  may  write 


%y)  =  I 


^{y-yi)  ifxo^x^xi 
0  otherwise. 


The  CRT  of  the  horizontal  line  segment  is  then  given  by  the  expression 

if(w-a;i)2  +  y?  ^  (n-a;o)2  +  y? 


if  u  ^  Xi: 


if  u  ^  xq: 


/  R  ' 

1  0  c 

/  R  ' 

1  0  c 


otherwise 


if  (u  -  xq)^  +  yf  ^{u-  xi)^  +  y? 
otherwise 


9{u,t)  = 


0  if  t  ^  yi 

2r  +  yf  ^  +  yf  > 

ifxo<u<xi:  if  (u  -  xo)^  +  yf  <  and  (u  -  xi)^  +  y?  ^ 

1^  if  (u  -  Xo)^  +  yf  >  f  and  {u  -  xif  +  yf  <  f 

0  if  (u  -  xo)^  +  y‘i  <  and  (u  -  xi)^  +  yf  < 


using  the  density  at  the  intersection  calculated  in  Section  7.2.6.  Figure  11  shows  the 
geometric  configurations  that  contribute  to  the  cases  in  (75).  Figure  12  shows  a  density 
plot  of  the  CRT  of  the  horizontal  line  segment. 
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0) 


(k) 


(1) 


Figure  11:  The  various  cases  under  which  a  circle  can  intersect  a  horizontal  line  segment. 


7.2.8  Simple  Polynomials  and  Other  Functions 


Table  1  presents  the  CRT  and  SORT  of  a  number  of  simple  target  functions.  Note 
that  Si{-)  is  the  sine  integral,  defined  to  be 


dt. 


7.3  Analytic  Inverse  Transforms 

In  this  section  we  will  use  the  resampling  transforms  presented  in  Section  6.1  to  dc'rive 
the  inverse  CRT  for  target  functions  with  known  CRT. 
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Figure  12:  The  CRT  of  a  horizontal  line  segment  yi  =  2  in  the  interval  xq  =  — 2,xi  =  2. 
A  square-root  mapping  to  compress  the  grey  scale  has  been  used  in  this  density  plot  to 
make  the  “arms”  visible. 


7.3.1  Bessel  Function 

Let  our  target  function  be  the  CRT  of  the  Gaussian  radial  basis  function  in  (72) , 

g{u,  t)  =  27rte-’^(*'+“'>/o(27rtu).  (76) 

Firstly,  note  that  g{u,t)  is  odd  in  the  second  variable  t,  so  we  assume  that  g{u,t)  =  0  for 
t  <  0.  Then  following  the  steps  presented  in  Section  6.1,  from  (52)  and  (53), 

/OO 

•OO 

fOO 

=  27re-’^“  /  te^^  ^^^-'^'>M2mtu)dt, 

Jo 


then  using  (6.631.4)  from  [13], 


Then  from  (54), 


G{u,u) 


2 

27rtu>u 
^  1— 2ia; 


1  —  2iu} " 


«(<:■.«  =  e' 


e  l-2i0 

l-2iP' 
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f{x,y) 

(7^c/)(u,  t) 

a 

lant 

ant 

ax 

2a'Ktu 

antu 

ay 

0 

2at^ 

ax"^ 

a-Ktit^  +  2u^) 

\ant{t^  +  2n^) 

ay"^ 

arct^ 

\ant^ 

ax^  +  by"^ 

Tc[at^  +  bt^  +  2atu^) 

\n{at^  +  bt^  +  2atu^) 

axy 

0 

2at^u 

ax^ 

a'rrt{St^u  +  2u^) 

\ant{Zt^u  +  2u^) 

ay^ 

0 

lat^ 

ax^y 

airt^u 

^ant^u 

axy"^ 

0 

at^{^t'^  +  2u^) 

a 

X 

2aiTt  /  u+t 
t+u  Y  u—t 

ant  f  u+t 
t+u  Y  u—t 

a 

2aiTtu  /  u—t 

airtu  u—t 

sin(7r^-f ) 
n9-^ 

lt—u)'^{t+u)  Y  'u.+t 

(i— u)^(i+u)  Y  u+t 

Table  1:  The  CRTs  and  SCRTs  of  a  number  of  simple  functions. 


Then  from  (55), 


k{x,r) 


,0)  da  d^ 


g— 27ri;Sr 
g-27ri/9r 

l-2iP 


— 

-27riax  ^  ^ 

1-2*^ 

g-7r(2iax+jfj3) 


d^, 
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then  using  (3.323.2)  from  [13], 


Using  (56), 


/oc  . 

00  ^—2wi0r—nx'^{l—2iP) 


-L 


oo  v/l  -  2ip 

oo  p— 27ri^(r— x^) 


dp 


i 


-2  e' 


00  \/l  -  2ip 

7r(i — x^) 


=  e 


r  — 


k{x,r)  — 


\/t  —  x^ 


f{x,y)  =  yk{x,x^  +  y^ 


which  agrees  with  (72). 


8  Options  for  Inverting  the  Circular  Radon 

Transform 

In  this  report  so  far  we  have  presented  two  equivalent  methods  in  Section  6  for  inverting 
the  CRT  via  resampled  Fourier  transforms.  In  this  section  we  discuss  a  number  of  options 
by  which  this  inversion  could  be  implemented.  In  addition,  we  present  ideas  for  another 
novel  method  that  uses  an  approximation  to  the  fast-time  matched  filtered  radar  signal 
via  a  linear  combination  of  analytic  circular  Radon  transformed  signals. 


8.1  Resampled  Fourier  Methods 

In  Section  6.1,  we  presented  a  succession  of  resamplings  and  Fourier  transforms  that 
invert  the  CRT.  There  are  a  number  of  different  ways  of  implementing  this  inversion, 
which  we  will  now  discuss. 

1.  The  most  obvious  approach  to  take  is  to  undertake  the  implementations  using  the 
FFT  algorithm  in  conjunction  with  interpolation  to  a  uniform  grid  after  each  re¬ 
sampling.  This  is  a  relatively  simple  approach,  but  not  particularly  efficient  because 
interpolation  is  generally  computationally  intensive. 

2.  Two  of  the  three  interpolation  steps  could  be  eliminated  by  using  non-uniformly 
sampled  FFT  algorithms  [3,  4,  11,  29,  33,  37],  which  have  computational  complexity 
of  the  order  of  the  standard  FFT. 
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3.  A  combination  of  linear  approximation  and  analytic  expressions  can  be  used  to 
remove  the  need  to  compute  some  of  the  resampled  Fourier  transforms.  For  example, 
in  the  first  step  of  Section  6.1,  we  compute 

/OO 

-CX) 

If  we  approximate  g{u,  t)  by  a  sum  of  Gaussian  radial  basis  functions. 


9{u,t) 


E 


Cj€.  '  2^ 


where  for  the  moment  we  assume  that  cr^  =  1  /2,  then  we  can  write 


G{u,u)  j 


-n{{t -ti)^ +{u-Ui)'^)  27Tiujt'^ 


dt 


n 


=  E 


cr 


v/l  —  2iuj 


Following  the  next  step 

'■OO  /•OO 


/OO  roo 
■OO  J  “OO 

=/7 

-n 

J  — OO  J  - 


/?)  da  d/3 

L\J 


-2T:icxx 


L 


—  OO  J  -OO 
OO 


E' 


-  2i[3 


7rf„2  +  iii£  .  iiL^L  4.  \ 


da  dp 


>  c 


/Tpr 

i=i  y^  +  iF 


2^l(p-p  +  2’ii^)0  +  H2p+2uJ  -Au,x+2i'^)0'^+4t'^0^ 

:€  1+4/?  2 


dp. 

(77) 


where  the  integral  (77)  cannot  be  computed  analytically  and  so  must  be  carried  out 
by  some  numerical  means  such  £is  quadrature  or  by  sampling  the  Fourier  transform 
integrand  and  using  the  FFT.  There  are  many  other  possible  choices  of  basis  function 
apart  from  the  Gaussian  radial  basis  function  used  in  this  example,  and  many  other 
combinations  of  analytical  and  numerical  steps  could  be  employed  to  implement  the 
resampled  Fourier  methods  of  Section  6. 


8.2  Linear  Combination  Method 

In  this  section,  we  develop  an  idea  for  SAR  image  formation  that  is  tlie  logical  extension 
of  the  third  idea  presented  in  the  previous  subsection.  The  idea  is  to  ai)proxiinate  the 
fast-time  slow-time  signal  as  a  linear  combination  of  the  transformed  basis,  and  then  the 
reconstructed  image  is  the  same  weighted  sum  of  the  original  basis  because  the  CRT  is  a 
linear  transform.  The  practical  impediments  to  this  approach  are  that  in  a  real  system, 
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only  a  portion  of  the  transformed  signal  is  observed  because  of  the  finite  extent  of  the 
antenna  beam  pattern.  The  portion  of  transformed  basis  function  observed  will  depend 
on  the  angle  to  the  centre  of  the  swath  and  the  antenna  pattern. 

We  start  out  by  considering  the  Gaussian  RBF  transform  pair, 


f{x,y)  =  e  2-" 


9iu,t)  =  {'R,cf){u,i) 


p27r 

/ 

Jo 


u  +  t  cos  9,  t  sin  6)t  d9 


i2+„2+j,2  /  tJv?  +yf 


(78) 


where  a  is  the  width  of  the  Gaussian  RBF.  Now  (78)  describes  an  hyperbolic  shape  with 
a  finite  cross-section  (an  hyperbolic  ensemble  (Figure  9)),  which  would  approach  (70),  the 
CRT  of  a  delta  function,  in  the  limit  as  a  goes  to  zero  (to  a  constant  multiple).  We  need  to 
determine  the  portion  of  the  hyperbolic  ensemble  that  is  relevant  to  a  particular  imaging 
geometry  and  antenna  combination. 

The  average  squint  angle  (off  trajectory)  of  a  target  and  its  observed  angular  extent 
(dictated  by  the  distance  to  target,  beam  width,  and  antenna  steering)  will  determine  the 
interval  [uo,ui]  over  which  the  CRT  of  the  target  will  be  observed.  If  the  target  area  is 
small  enough,  this  can  be  assumed  to  be  constant  for  all  targets  across  the  target  area. 
Next,  we  determine  the  weighted  sum  of  translates  of  the  windowed  version  of  (78)  that 
best  matched  the  received  signal.  Finally,  the  reconstructed  image  is  given  by  the  same 
weighted  sum  of  the  same  translates  of  the  Gaussian  RBF.  There  are  many  details  yet  to 
be  sorted  out,  however,  to  make  this  approach  a  practical  one. 

Note  that  in  this  suggested  approach,  we  have  in  broad  terms  substituted  the  stationary 
phase  approximation  of  the  RMA  by  another.  However,  it  would  be  possible  to  control  the 
fidelity  of  the  new  approximation  to  achieve  any  desired  degree  of  accuracy  in  a  natural 
way,  something  that  is  not  possible  in  the  stationary  phase  approximation. 


9  Motion  Compensation 

We  will  now  examine  how  motion  compensation  can  be  incorporated  into  the  schemes 
presented  here  for  image  formation  via  inversion  of  the  CRT. 

Let  us  examine  the  effect  of  perturbing  a  sample  point  on  the  transform  space  data. 
This  is  one  step  towards  our  goal  of  developing  an  inversion  process  for  arbitrary  imaging 
trajectory.  Figure  13  shows  the  geometry  of  the  problem.  We  consider  displacement 
perpendicular  to  the  imaging  trajectory  only,  because  displacement  along  the  imaging 
trajectory  is  simple  to  deal  with  as  we  discuss  in  a  moment  (but  also  see  Section  7.1). 
From  this  diagram  we  can  see  the  effect  of  the  perturbation,  which  is  to  change  the  radius 
of  the  intersecting  circle  at  the  sample  point  x  =  u  along  the  imaging  axis. 

We  can  calculate  the  effect  upon  the  intersecting  circle.  Let  d  denote  the  displacement 
from  the  imaging  axis  as  in  Figure  13.  Then  the  error  in  t  of  the  transform  of  the  points 
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Figure  IS:  The  geometry  of  a  displaced  point. 


in  the  plane  is  given  by 

t  —  t'  =  \/{x  -  uY  +  2/2  -  \J{x  -  +  {y  -  df^ 

—  t  -  \/d?‘  +  —  2dt  sin  0.  (79) 


If  we  displace  the  points  in  the  target  plane  along  the  direction  9  to  t',  then  the  grid  has 
been  distorted  to  give  that  shown  in  Figure  14. 


Motion  compensation  for  reconstruction  via  backprojection  is  simply  a  matter  of  ad¬ 
justing  the  coordinates  to  compensate  for  the  platform  movement, 


2^^{Xn  -U-  Xeju))^  -h  {pn  - 

"  ) 


du. 


(80) 


More  generally,  motion  compensation  can  be  undertaken  as  follows  [32].  Let  our  sensor 
platform  trajectory  be  given  by  [u  -I-  Xe(u),yeiu)].  Note  that  in  standard  approaches  to 
image  formation,  standard  FFT  routines  are  used  that  require  uniform  sampling.  Conse¬ 
quently,  motion  compensation  in  these  cases  must  account  for  an  error  in  the  along-track 
direction  to  compensate  for  platform  motion  causing  deviation  from  uniform  slow-time 
sampling.  However,  if  non-uniformly  sampled  FFT  routines  are  used  then  this  sampling  is 
not  a  problem  so  motion  compensation  can  be  treated  as  an  error  in  only  the  cro.ss-track 
coordinate  y. 


Let  our  generic  SAR  signal  be  given  by 


s(u 


,t)  =  an  p  i  - 


2y/{Xn  -  u)2  +  yl 


n  \  / 

Then  using  the  shift  theorem,  the  Fourier  transform  of  the  generic  SAR  signal  with  respect 
to  fast- time  t  is 


S{u.u>)  =P(w)^a„  exp  (—^■Kik\/{xn  -  uY  +  y 


where  k  =  uifc  and  is  called  the  wavenumber.  Consequently,  the  measured  SAR,  signal  of 
target  {xn,yn)  di  the  (u.o;)  domain  under  motion  errors  is 

5„(n,w)  =  exp  ^-47riA,V (a^,  -  u  -  Xc.(u)f  +  {y„  -  y,.(u))~'^ 
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Figure  14:  The  distortion  of  the  cartesian  grid  due  to  the  displaced  point. 


where  the  total  measured  SAR  signal  in  the  (u,  uj)  domain  is 

S{U,U)  =  P{u)'Y^anSn{u,U}). 

n 

We  can  rewrite  this  as 

5n(u,w)  =  aen{xn  “  u,yn,oj)exp  (^-iniky/ (Xn  -  u)^  + 
where  the  motion  phase  error  function  is  given  by 

aen{xn  “  =  exp(47ri/ere„(u)) 

and  the  radial  error  for  the  n-th  target  is 

ren{u)  =  -y/{Xn  -U-  Xe{u)Y  +  (y„  -  ye(^i))^  +  {Xn  “  u)‘^  +  yl- 


Now,  provided  the  fluctuations  of  aen{xn  —  u,yni<^)  are  small  compared  with  the  fluc¬ 
tuations  of  exp  ^—4niky/{xn  —  u)‘^  -I-  y^ ,  then  by  the  method  of  stationary  phase  ([32], 
p.  99)  the  slow-time  Fourier  transform  of  5„(n,a;)  is 

Sen{ku,<jo)  =  Aeniku,oj)exp  (^-2nikuXn  -  2'Ki\J\k'^  -  kl  ynj  , 

where 


Agri  (2/j  cos  ^72  (n) ,  Cc^)  —  a^i^i^Xji  u.^yyi^Lo') 
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and  the  aspect  angle  of  the  n-th  target  when  the  radar  is  located  at  {u,  0)  is  given  by 

6n{u)  =  tan"^  ^  . 

\Xn  -  uj 

Now  from  the  SAR  mapping 

kyiku^u)  =  v/4A:2  -  kl 

and  the  slow-time  to  slow-time  frequency  domain  mapping 

ku  =  2k  cos  6n{u) 

we  obtain  that 

2k=^kl  +  kl 

U  =  Xn-  ^Vn- 

Ky 

Therefore,  ~  u,  directly  maps  into  the  {ky^,u})  domain,  and  furthermore,  the 

{ku,u))  domain  directly  maps  into  the  target  spatial  frequency  domain  {kx,ky).  Conse¬ 
quently,  provided  that  the  fluctuations  of  aen{xn  -  u,yn-,ijj)  are  small  compared  with  the 
motion  error- free  signal,  the  motion  phase  error  function  —  u,  can  be  modelled 

as  a  filter  in  the  spatial  frequency  domain  {kx,ky).  This  filter,  denoted  Hen{kx,  ky),  varies 
spatially  with  the  coordinates  of  the  target,  and  is  given  by 

Hen{kxiky)  —  QgfilXn  —  U,yii,Uj) 

=  exp  {Airik  re„(u)) .  (81) 


For  narrow  beamwidth  SAR  systems,  the  radial  motion  error  function  is  approximated 

by 

feniu)  ~  Xe{u)  COS  6c  +  yeiu)  Sm  9c 

where  dc  is  the  average  squint  angle  (off  trajectory)  of  the  target  area.  Therefore,  narrow 
beamwidth  motion  compensation  is  performed  via 

Sn{u,  u)  exp  (—47rikxe(u)  cos  6c  —  Anikydu)  sin  0c) ,  (82) 


which  in  the  broadside  case  {6c  =  7r/2)  simplifies  to 

-\/{xn-u-  Xe{u))'^  +  (y„  -  a;e(u))2  4-  yj{xn  -  n)^  +  ylSn{u,u)eyi\>{-Amkyc{u)) , 
a  filter  that  is  independent  of  the  position  of  the  coordinates  of  the  target. 

In  wide  beamwidth  systems,  the  assumption  that  the  fluctuations  of  aen{xn  -  u,  Un^ix) 
are  small  compared  with  that  of  the  motion  error-free  signal  is  violated.  One  sohition  is 
to  perform  narrow  beamwidth  compensation  first  using  (82),  which  reduces  the  dynamic 
range  of  the  errors,  and  then  apply  the  shift- varying  filter  (81)  to  compensate  for  the 
remaining  positional  errors.  After  the  narrow  beam  compensation  to  average  squint  angle 
6c,  the  motion  error  still  to  be  compensated  for  is  given  by  the  expression 

Vexy  =  -Vixn  -ti-  Xeiu))'^  +  (?y„  -  Xe{u))^  +  'J (x,,  -  u)^  -f  _  Xc{u)  COS  6c  -  yc{u)sm6c. 
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10  Conclusion 

The  purpose  of  this  report  has  been  to  present  an  alternative  approach  to  SAR  image 
formation  from  the  perspective  of  circular  Radon  transforms  (CRTs).  This  approach  can 
naturally  handle  high-squint  and  ultra-wideband  SAR.  We  have  reviewed  the  literature 
on  the  inversion  of  circular  Radon  transforms,  some  of  which  has  been  undertaken  by 
SAR  researchers  in  Sweden  for  the  CARABAS-II  sensor.  We  present  a  number  of  novel 
possible  approaches  to  the  CRT  inversion  problem  that  could  be  developed  into  practical 
SAR  image  formation  algorithms,  and  we  will  undertake  this  development  in  follow-on 
research. 
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